home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 20 / Cream of the Crop 20 (Terry Blount) (1996).iso / os2 / xdsn217.zip / SAMPLES / SIMPLE / linnew.c < prev    next >
C/C++ Source or Header  |  1996-06-04  |  23KB  |  881 lines
  1. /*
  2. **
  3. ** LINPACK.C        Linpack benchmark, calculates FLOPS.
  4. **                  (FLoating Point Operations Per Second)
  5. **
  6. ** Translated to C by Bonnie Toy 5/88
  7. **
  8. ** Modified by Will Menninger, 10/93, with these features:
  9. **  (modified on 2/25/94  to fix a problem with daxpy  for
  10. **   unequal increments or equal increments not equal to 1.
  11. **     Jack Dongarra)
  12. **
  13. ** - Defaults to double precision.
  14. ** - Averages ROLLed and UNROLLed performance.
  15. ** - User selectable array sizes.
  16. ** - Automatically does enough repetitions to take at least 50 CPU seconds.
  17. ** - Prints machine precision.
  18. ** - ANSI prototyping.
  19. **
  20. ** To compile:  cc -O -o linpack linpack.c -lm
  21. **
  22. **
  23. */
  24.  
  25. #include <stdio.h>
  26. #include <stdlib.h>
  27. #include <math.h>
  28. #include <time.h>
  29. #include <float.h>
  30.  
  31. #define DP
  32.  
  33. #ifdef SP
  34. #define ZERO        0.0
  35. #define ONE         1.0
  36. #define PREC        "Single"
  37. #define BASE10DIG   FLT_DIG
  38.  
  39. typedef float   REAL;
  40. #endif
  41.  
  42. #ifdef DP
  43. #define ZERO        0.0e0
  44. #define ONE         1.0e0
  45. #define PREC        "Double"
  46. #define BASE10DIG   DBL_DIG
  47.  
  48. typedef double  REAL;
  49. #endif
  50.  
  51. static REAL linpack  (long nreps,int arsize);
  52. static void matgen   (REAL *a,int lda,int n,REAL *b,REAL *norma);
  53. static void dgefa    (REAL *a,int lda,int n,int *ipvt,int *info,int roll);
  54. static void dgesl    (REAL *a,int lda,int n,int *ipvt,REAL *b,int job,int roll);
  55. static void daxpy_r  (int n,REAL da,REAL *dx,int incx,REAL *dy,int incy);
  56. static REAL ddot_r   (int n,REAL *dx,int incx,REAL *dy,int incy);
  57. static void dscal_r  (int n,REAL da,REAL *dx,int incx);
  58. static void daxpy_ur (int n,REAL da,REAL *dx,int incx,REAL *dy,int incy);
  59. static REAL ddot_ur  (int n,REAL *dx,int incx,REAL *dy,int incy);
  60. static void dscal_ur (int n,REAL da,REAL *dx,int incx);
  61. static int  idamax   (int n,REAL *dx,int incx);
  62. static REAL second   (void);
  63.  
  64. static void *mempool;
  65.  
  66.  
  67. void main(void)
  68.  
  69.     {
  70.     char    buf[80];
  71.     int     arsize;
  72.     long    arsize2d,memreq,nreps;
  73.     size_t  malloc_arg;
  74.  
  75.     while (1)
  76.         {
  77.         printf("Enter array size (q to quit) [200]:  ");
  78.         fgets(buf,79,stdin);
  79.         if (buf[0]=='q' || buf[0]=='Q')
  80.             break;
  81.         if (buf[0]=='\0' || buf[0]=='\n')
  82.             arsize=200;
  83.         else
  84.             arsize=atoi(buf);
  85.         arsize/=2;
  86.         arsize*=2;
  87.         if (arsize<10)
  88.             {
  89.             printf("Too small.\n");
  90.             continue;
  91.             }
  92.         arsize2d = (long)arsize*(long)arsize;
  93.         memreq=arsize2d*sizeof(REAL)+(long)arsize*sizeof(REAL)+(long)arsize*sizeof(int);
  94.         printf("Memory required:  %ldK.\n",(memreq+512L)>>10);
  95.         malloc_arg=(size_t)memreq;
  96.         if (malloc_arg!=memreq || (mempool=malloc(malloc_arg))==NULL)
  97.             {
  98.             printf("Not enough memory available for given array size.\n\n");
  99.             continue;
  100.             }
  101.         printf("\n\nLINPACK benchmark, %s precision.\n",PREC);
  102.         printf("Machine precision:  %d digits.\n",BASE10DIG);
  103.         printf("Array size %d X %d.\n",arsize,arsize);
  104.         printf("Average rolled and unrolled performance:\n\n");
  105.         printf("    Reps Time(s) DGEFA   DGESL  OVERHEAD    KFLOPS\n");
  106.         printf("----------------------------------------------------\n");
  107.         nreps=1;
  108.         while (linpack(nreps,arsize)<40.)
  109.             nreps*=2;
  110.         free(mempool);
  111.         printf("\n");
  112.         }
  113.     }
  114.  
  115.  
  116. static REAL linpack(long nreps,int arsize)
  117.  
  118.     {
  119.     REAL  *a,*b;
  120.     REAL   norma,t1,kflops,tdgesl,tdgefa,totalt,toverhead,ops;
  121.     int   *ipvt,n,info,lda;
  122.     long   i,arsize2d;
  123.  
  124.     lda = arsize;
  125.     n = arsize/2;
  126.     arsize2d = (long)arsize*(long)arsize;
  127.     ops=((2.0*n*n*n)/3.0+2.0*n*n);
  128.     a=(REAL *)mempool;
  129.     b=a+arsize2d;
  130.     ipvt=(int *)&b[arsize];
  131.     tdgesl=0;
  132.     tdgefa=0;
  133.     totalt=second();
  134.     for (i=0;i<nreps;i++)
  135.         {
  136.         matgen(a,lda,n,b,&norma);
  137.         t1 = second();
  138.         dgefa(a,lda,n,ipvt,&info,1);
  139.         tdgefa += second()-t1;
  140.         t1 = second();
  141.         dgesl(a,lda,n,ipvt,b,0,1);
  142.         tdgesl += second()-t1;
  143.         }
  144.     for (i=0;i<nreps;i++)
  145.         {
  146.         matgen(a,lda,n,b,&norma);
  147.         t1 = second();
  148.         dgefa(a,lda,n,ipvt,&info,0);
  149.         tdgefa += second()-t1;
  150.         t1 = second();
  151.         dgesl(a,lda,n,ipvt,b,0,0);
  152.         tdgesl += second()-t1;
  153.         }
  154.     totalt=second()-totalt;
  155.     if (totalt<0.5 || tdgefa+tdgesl<0.2)
  156.         return(0.);
  157.     kflops=2.*nreps*ops/(1000.*(tdgefa+tdgesl));
  158.     toverhead=totalt-tdgefa-tdgesl;
  159.     if (tdgefa<0.)
  160.         tdgefa=0.;
  161.     if (tdgesl<0.)
  162.         tdgesl=0.;
  163.     if (toverhead<0.)
  164.         toverhead=0.;
  165.     printf("%8ld %6.2f %6.2f%% %6.2f%% %6.2f%%  %9.3f\n",
  166.             nreps,totalt,100.*tdgefa/totalt,
  167.             100.*tdgesl/totalt,100.*toverhead/totalt,
  168.             kflops);
  169.     return(totalt);
  170.     }
  171.  
  172.  
  173. /*
  174. ** For matgen,
  175. ** We would like to declare a[][lda], but c does not allow it.  In this
  176. ** function, references to a[i][j] are written a[lda*i+j].
  177. */
  178. static void matgen(REAL *a,int lda,int n,REAL *b,REAL *norma)
  179.  
  180.     {
  181.     int init,i,j;
  182.  
  183.     init = 1325;
  184.     *norma = 0.0;
  185.     for (j = 0; j < n; j++)
  186.         for (i = 0; i < n; i++)
  187.             {
  188.             init = (int)((long)3125*(long)init % 65536L);
  189.             a[lda*j+i] = (init - 32768.0)/16384.0;
  190.             *norma = (a[lda*j+i] > *norma) ? a[lda*j+i] : *norma;
  191.             }
  192.     for (i = 0; i < n; i++)
  193.         b[i] = 0.0;
  194.     for (j = 0; j < n; j++)
  195.         for (i = 0; i < n; i++)
  196.             b[i] = b[i] + a[lda*j+i];
  197.     }
  198.  
  199.  
  200. /*
  201. **
  202. ** DGEFA benchmark
  203. **
  204. ** We would like to declare a[][lda], but c does not allow it.  In this
  205. ** function, references to a[i][j] are written a[lda*i+j].
  206. **
  207. **   dgefa factors a double precision matrix by gaussian elimination.
  208. **
  209. **   dgefa is usually called by dgeco, but it can be called
  210. **   directly with a saving in time if  rcond  is not needed.
  211. **   (time for dgeco) = (1 + 9/n)*(time for dgefa) .
  212. **
  213. **   on entry
  214. **
  215. **      a       REAL precision[n][lda]
  216. **              the matrix to be factored.
  217. **
  218. **      lda     integer
  219. **              the leading dimension of the array  a .
  220. **
  221. **      n       integer
  222. **              the order of the matrix  a .
  223. **
  224. **   on return
  225. **
  226. **      a       an upper triangular matrix and the multipliers
  227. **              which were used to obtain it.
  228. **              the factorization can be written  a = l*u  where
  229. **              l  is a product of permutation and unit lower
  230. **              triangular matrices and  u  is upper triangular.
  231. **
  232. **      ipvt    integer[n]
  233. **              an integer vector of pivot indices.
  234. **
  235. **      info    integer
  236. **              = 0  normal value.
  237. **              = k  if  u[k][k] .eq. 0.0 .  this is not an error
  238. **                   condition for this subroutine, but it does
  239. **                   indicate that dgesl or dgedi will divide by zero
  240. **                   if called.  use  rcond  in dgeco for a reliable
  241. **                   indication of singularity.
  242. **
  243. **   linpack. this version dated 08/14/78 .
  244. **   cleve moler, university of New Mexico, argonne national lab.
  245. **
  246. **   functions
  247. **
  248. **   blas daxpy,dscal,idamax
  249. **
  250. */
  251. static void dgefa(REAL *a,int lda,int n,int *ipvt,int *info,int roll)
  252.  
  253.     {
  254.     REAL t;
  255.     int idamax(),j,k,kp1,l,nm1;
  256.  
  257.     /* gaussian elimination with partial pivoting */
  258.  
  259.     if (roll)
  260.         {
  261.         *info = 0;
  262.         nm1 = n - 1;
  263.         if (nm1 >=  0)
  264.             for (k = 0; k < nm1; k++)
  265.                 {
  266.                 kp1 = k + 1;
  267.  
  268.                 /* find l = pivot index */
  269.  
  270.                 l = idamax(n-k,&a[lda*k+k],1) + k;
  271.                 ipvt[k] = l;
  272.  
  273.                 /* zero pivot implies this column already
  274.                    triangularized */
  275.  
  276.                 if (a[lda*k+l] != ZERO)
  277.                     {
  278.  
  279.                     /* interchange if necessary */
  280.  
  281.                     if (l != k)
  282.                         {
  283.                         t = a[lda*k+l];
  284.                         a[lda